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Recent theoretical advances have established that the electric polarization in an insulating crystal 
can be viewed as a multivalued quantity that is determined by certain Berry phases associated with 
the occupied Bloch bands. The application of this approach to the computation of piezoelectric coef- 
ficients is not entirely straightforward, since a naive determination of the ( "improper" ) piezoelectric 
coefficients from finite differences of the polarization at nearby strain states leads to a dependence 
upon the choice of "branch" of the polarization. The purpose of the present paper is to clarify that 
if one calculates instead the "proper" piezoelectric response, the branch dependence is eliminated. 
From this analysis, a simplified recipe for the direct finite-difference computation of the proper 
piezoelectric coefficients emerges naturally. 
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I. INTRODUCTION 

The calculation of spontaneous polarization and piezo- 
electric response within the framework of first-principles 
methods of electronic structure theory has proven to be 
a rather subtle problem. In a landmark paper, Martin 
jl| showed that the piezoelectric tensor is well-defined as 
a bulk quantity in a crystalline insulator. However, at 
that time it was far from clear whether the spontaneous 
polarization itself could be regarded as a bulk property 
in the same sense, and calculations of piezoelectric con- 
stants by finite differences of spontaneous polarization 
were therefore not possible. 

The situation changed in 1993 with the development 
of the "Berry-phase" theory of polarization [0J|] , which 
provided a direct and straightforward method for com- 
puting the electric polarization. (For a useful review, see 
Ref. Q].) Nevertheless, some subtleties remain regarding 
the computation of the piezoelectric tensor components 
by finite differences ||[|. First, the Berry-phase theory 
gives the polarization as a multivalued quantity, and the 
piezoelectric response that would be computed from a 
given one of the many branches is not invariant with re- 
spect to choice of branch. Second, a distinction is made 
between the "proper" and "improper" piezoelectric re- 
sponse , and it might not be clear which of these is 
to be associated with the finite-difference calculation. 

The purpose of the present paper is to elucidate the 
physics of the spontaneous polarization, the piezoelec- 
tric response, and the relations between the two. It is 
clarified that the improper polarization is the one given 
by the naive finite-difference approach, and that while 
this quantity is indeed branch-dependent, the proper po- 
larization, which should be compared with experiment, 
is not. As a result of this analysis, a simplified recipe 
for the direct finite-difference computation of the proper 
piezoelectric response is given. 



II. BERRY-PHASE THEORY OF POLARIZATION 

We consider a periodic insulating crystal in zero macro- 
scopic electric field, and assume that the electronic 
ground state can be described by a one-electron Hamilto- 
nian H as in density-functional or Hartree-Fock theory. 
The eigenstates of H are the Bloch functions ijj n k with 
energies e„k, and it is conventional to define the cell- 
periodic Bloch functions 



Mnk(r) = e lk ' r V«k(r) 



(1) 



having periodicity u n k( r ) = w«k(R- + r), where R is any 
lattice vector. The contribution of the n'th occupied 
band to the spontaneous electric polarization of the crys- 
tal can then be written 10.01 
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d 3 k (u„k|V k |u Ilk ) 



(2) 



We take the convention that n runs over bands and spin, 
so a factor of two would need to be inserted in Eq. (||) 
to account for paired spins. The total spontaneous po- 
larization is then given by 



it 



(3) 



where Z T and r T are the atomic number and cell position 
of the r'th nucleus in the unit cell, and fl is the unit cell 
volume. 

Strictly speaking, Eq. @ applies only to an isolated 
band, i.e., a band for which e ra k does not become de- 
generate with any other band at any point in the Bril- 
louin zone. This restriction is not essential; methods for 
extending the analysis to composite groups of occupied 
bands containing arbitrary degeneracies and crossings 
have been developed as described in Refs. [§-[f). How- 
ever, for simplicity of presentation, it will be assumed 
here that only isolated bands are present. For the same 
reason, spin degeneracy is suppressed throughout. 

There is a certain arbitrariness inherent in Eq. (^) as- 
sociated with the freedom to choose the phases of the 



1 



Bloch functions f/W- For, suppose we make a different 
choice 

|^nk) = e^ 0t) |^nk) • (4) 

We shall refer to this as a "gauge transformation" of 
the Bloch functions. Note that the choice of /3(k) is re- 
stricted by the fact that k and k + G label the same 
wavefunction (where G is a reciprocal lattice vector), so 
that f3(k + G) — /3(k) must be an integer multiple of 2ir 
for any G. Thus, the most general form of /3(k) is 

/?(k)=/? per (k)+k-R (5) 

where /3 por is a periodic function in fc-space and R is 
some real-space lattice vector. Letting P„ be the result 
of inserting the w„k in place of the u n k in Eq. (||), and 
using Eqs. (Q), (f|), and (||), one finds 

P„ =Pn-f • (6) 

Thus, while the contribution of this band to the electronic 
polarization is not absolutely gauge-invariant, it is gauge- 
invariant modulo e/n times a real-space lattice vector. 
Actually, this is precisely the type of qualified invariance 
we should have expected. After all, the choice of the 
location r T of the atom representing sublattice r in the 
unit cell has a similar ambiguity; we could just as well 
choose r r = r r + R', where R' is another lattice vector, 
leading to precisely the same kind of "modulo eR/fi" 
ambiguity in the expression for P in Eq. (||). 

Perhaps the most natural way to incorporate this kind 
of ambiguity in the definition of the polarization is to 
regard P as a multivalued quantity; that is, it simulta- 
neously takes on a lattice of values given by some P^ 
(here is a "branch" label) and all its periodic images 
p( 6 ' + eR/fi (with R running over all lattice vectors of 
the crystal). To interpret this intuitively, we can say 
that from the point of view of its dipolar properties, the 
real insulator behaves like a fictitious crystal composed 
of two sublattices of point ±e charges, with the sublat- 
tice of — e charges displaced relative to the +e sublattice 
by -OP/e. That is, choosing one of the +e charges as 
the origin, — fiP/e takes on a lattice of values that is 
precisely the lattice of positions of the — e charges. 

This is illustrated in Fig. 1 for an imaginary tetragonal 
crystal (dimensions a x a x c) with one monovalent ion 
located at the cell corners, and a single (spinless) elec- 
tron band giving rise to the distributed electron charge 
indicated schematically by the contours in Fig. 1(a). (We 
assume that M z mirror symmetry is broken in some way.) 
Eq. (|J) then gives the location 

r n = -OP n /e (7) 

of the effective unit point charge — e illustrated in 
Fig. 1(b). As discussed in Refs. ^-j4|, this location is 




FIG. 1. Illustrative tetragonal crystal (cell dimensions 
a x a x c) having one monovalent ion at the cell corner (ori- 
gin) and one occupied valence band, (a) The distributed 
quantum-mechanical charge distribution associated with the 
electron band, represented as a contour plot, (b) The dis- 
tributed electron distribution has been replaced by a unit 
point charge — e located at the Wannier center r n , as given 
by the Berry-phase theory. 

just the charge center of the Wannier function associated 
with the electron band. The polarization will then take 
on a lattice of values having x, y, and z components of 
m\e/ac, rri2e/ac, and (7 + m^e/a 2 , respectively, where 
the n%i are integers. More generally, when several occu- 
pied bands are present, one can rewrite Eq. (||) as 

r n occ 

In practice, one proceeds by computing the component 
of P n along a particular crystallographic direction a via 
the quantity 

n 

0n.a G a ■ P n , (9) 

e 

where G Q is the primitive reciprocal lattice vector in di- 
rection a. In cases of simple symmetry (e.g., tetragonal 
or rhombohedral ferroelectric phases), a single <f> n suffices 
to determine P„, but in general P„ can be reconstructed 
from the three <^'s via 

P » = '^E^ R « - ( 10 ) 

a 

where R Q is the real-space primitive lattice vector cor- 
responding to G a . The 4> n ^ a are angle variables ("Berry 
phases") that are well-defined modulo 27r, given by 

4>n,a. = ^bz / d 3 fc (u„ k | - «G a • Vk|u„k) , (11) 
Jbz 

where 51bz = (27r) 3 /f2 is the Brillouin zone (BZ) volume. 

The 4> n ,a can be regarded as giving the position of the 
Wannier center for band n. For the toy crystal of Fig. 1, 
for example, and with the origin chosen at the cell corner, 
one would have <p x — <f> y — and <f> z = — 2^7. The 
practical calculation of the </> n ,a proceeds on a discrete 
mesh in reciprocal space, arranged as a two-dimensional 
grid of G a -oriented strings of k-points, as described in 
Refs. §f . 
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III. PIEZOELECTRIC RESPONSE 



The piezoelectric tensor of a crystal reflects the first- 
order change in spontaneous electric polarization in re- 
sponse to a first-order deformation of the crystal. The 
"improper" piezoelectric tensor is defined as |?],|| 



dPi 



c ljk 



dc 



in terms of the deformation 



drj = ^2 de jk rk 



(12) 



(13) 



where the symmetric and antisymmetric parts of dc rep- 
resent infinitesimal strains and rotations, respectively. 
On the other hand, the "proper" piezoelectric tensor can 
be defined as 



Cijk 



dJi 



dc 



(14) 



where J is the current density that flows through the bulk 
of the sample in adiabatic response to a slow deformation 
c = de/dt. According to the standard references the 
relation between the improper and proper piezoelectric 
tensors is 



Cijk = Cijk + SjkPi - SijPk 



(15) 



Writing out explicit tensor components, this last equa- 
tion becomes ft 



Czzz 


— Czzz 


i 


&ZXX 


— Czxx 


+ p z 


Czxy 


— Czxy 




C-zxz 


Czxz 




Czzx 


— C-zzx 


- p x 



(16) 



and similarly for permutations of the cartesian labels 
(but not for permutations of their position in the index 
triplet). It might seem strange at first sight that the ex- 
pressions for c zxz and c zzx have a different form, but this 
just reflects the fact that the deformation tensor e has 
been allowed to contain an antisymmetric part. 

Now in the Berry-phase theory, the polarization is a 
multivalued quantity, so that any particular value 
has to be identified by its branch label '6', and the cor- 
responding improper piezoelectric tensor is 



-ijk 



(fa) 



dc 



jk 



(17) 



Since P is well-defined modulo eR/f2, and both R and 
vary with the deformation e, Eq. ( |l7|) will clearly give 
different results for different choices of branch. This 



FIG. 2. Sketch of experiment for measuring proper piezo- 
electric coefficients. Strain e is applied to piezoelectric ma- 
terial (gray) sandwiched between grounded capacitor plates, 
and resulting current I is measured. 



branch-dependence is problematic; the piezoelectric ten- 
sor is measurable, and a suitable theory ought to give a 
unique value for it. 

Before proceeding, the reader is reminded that the 
piezoelectric response contains, in general, a "clamped- 
ion" part and an "internal-strain" part (5|^,^]. That is, 
one decomposes the actual deformation into a sum of 
two parts: a homogeneous strain in which the nuclear 
coordinates follow Eq. ( |l3| ) exactly (clamped- ion part), 
plus an internal distortion of the nuclear coordinates at 
fixed strain (internal-strain part). Since the latter oc- 
curs at fixed strain, all the subtleties about the branch- 
dependence and the proper-vs. -improper distinction dis- 
appear for this case. While the computation of the 
internal-strain part of the piezoelectric response may be 
tedious (requiring an iterative set of force calculations to 
determine the needed internal relaxations), it is straight- 
forward in principle. Consequently, for the remainder of 
this paper, the discussion refers to the clamped-ion re- 
sponse only unless explicitly stated otherwise. 



A. Branch-invariance of proper piezoelectric 
response 

While it is true that the improper piezoelectric re- 
sponse depends, in general, on choice of branch, it is in- 
stead the proper piezoelectric tensor that should be com- 
pared with experiment. Figure 2 shows a sketch of one 
possible experimental setup, in which a block of piezo- 
electric material is sandwiched between shorted conduct- 
ing electrodes, and the current / that flows in response to 
a deformation c is measured. As suggested by Eq. (|h!|), 
the proper piezoelectric response is related to the cur- 
rent that flows through the sample in response to the 
deformation, and is thus the experimentally measured 
quantity. Moreover, the induced current density j(r) is 
periodic with the lattice, so that its unit cell average J 
in Eq. ([u]) is perfectly well-defined, and consequently 
the proper piezoelectric tensor c cannot suffer from any 
dependence upon choice of branch. 

It is straightforward to check this branch-independence 
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of c explicitly. Since the polarizations for two different 
branch choices are related by 



p(6') = p(6) + £ R 



(18) 



one finds 



dPp = dP? ] - -Sr (KlRi + - dR, 



o: 2 

e 







so that 



(6') 



ijk ijk 

or, using Eq. (Ell), 







(19) 



(20) 



(*>') 



"ijk 



(6) 



(21) 



It is thus evident that Cijk as defined in Eq. ( |l5| ) is indeed 
independent of choice of branch. 

It is instructive to note that a similar argument applies 
to the part of the proper piezoelectric tensor arising from 
the ionic contribution P ion = (e/fi) ^2 T Z T r T in Eq. ||). 
Recalling that we are working in the clamped-ion ap- 
proximation, so that dr T follows the form of Eq. (13), 
one finds immediately that Ci on = by the same logic as 
for the previous paragraph. 

Indeed, the same logic would apply to Eq. (||) if the 
Wannier centers r„ would undergo a homogeneous defor- 
mation of the type ([k|) . In other words, the proper piezo- 
electric response is identically zero for a homogeneous 
deformation of both the ionic positions and the Wannier 
centers, in which case there is no charge flow through the 
interior of the crystal. 



B. Simplified finite-difference formula 

Of course, there is no reason to expect the Wannier 
centers r n to follow a homogeneous deformation, so c is 
not generally zero. But from this point of view, it be- 
comes evident that the the proper piezoelectric response 
is precisely a measure of the degree to which the Wan- 
nier centers fail to follow a homogeneous deformation. 
Or equivalently, returning to Eq. (|l0|), we see that the 
proper piezoelectric response measures just the variation 
of the Berry phases (f> n ,a with the strain deformation. 
More precisely, starting from Eqs. (|To|) , (|l2]), and (15), 
one finds 



Cijk 



"2tt n 



E 



dc 



R 



(22) 



Surface I 



+ + + + 



+ + + + 



+ + + + 



Surface II 



+ + + + 



+ + + + 



+ + + + 



FIG. 3. Two possible surface terminations of the lattice of 
point charges shown in Fig. 1(b). 



We have been working in the clamped-ion approximation, 
but in general, if there are internal relaxations accompa- 
nying the deformation, one can define a total Berry phase 
in direction a, 



so that 



7 - — - R 

rv J 



(23) 



(24) 



Naturally, the ionic contributions to d<j) a /de,jk vanish in 
the clamped-ion approximation. 

Equation (p2|), or its generalization (|2~i|), is the central 
result of this paper, and provides a simple and practical 
recipe for calculating the desired proper piezoelectric re- 
sponse. One simply computes the needed d<fi/ 1 de by finite 
differences, as (<// — </>)/(e' — e) for nearby strain config- 
urations e and e'. Then these d<j>/de are inserted into 
Eq. (22) or (E4I) to obtain the elements of the proper 



piezoelectric tensor. 



C. Relation to surface charges 



At the end of Sec. [II A, it was pointed out that a 



homogeneous deformation of the lattice of positive ionic 
and negative Wannier-center point charges would give 
rise to no internal current, and hence no proper piezo- 
electric response. This result can be made more intuitive 
by considering the connection between bulk polarization 
and surface charges ||. 

Consider, for example, a crystallite composed of N x 
N x N replicas of the unit cell shown in Fig. 1(b). In 
general there may be an arbitrariness in the choice of 
surface termination, as illustrated for the top surface of 
this crystallite in Fig. 3. For any given termination, the 
macroscopic surface charge density a is uniquely defined 
as J dz~p(z), where p~(z) is the average charge contained 
in a unit cell centered at vertical coordinate z (so that 
p vanishes either deep in the crystal or deep in the vac- 
uum and its integral is convergent). For the crystal of 



4 



Fig. 1, one finds a = je/a 2 and a = (7 — l)e/a 2 for the 
terminations of type I and II of Figs. 3(a) and 3(b), re- 
spectively. Referring back to Sec. ||, where it was found 
that P z = (7 + m^)e/a 2 , one confirms that the relation 
I 

a = P n (25) 

is satisfied for both terminations, the ambiguity of ter- 
mination corresponding to the choice of branch of P. 

For dcfinitcncss, assume that the surface terminations 
are such that the top and bottom surfaces of the crys- 
tallite have charge densities +je/a 2 and —je/a 2 on the 
top and bottom surfaces, respectively, and zero on the 
sides. Then the magnitude of the total charge on the 
top or bottom surface is just N 2 a 2 a = N 2 je, which is 
clearly independent of any homogeneous (7-preserving) 
deformation of the crystal. Thus, if this crystallite were 
inserted between grounded capacitor plates as in Fig. 2, 
no current would flow through the wire as a result of the 
homogeneous deformation. This is consistent with the 
vanishing of the proper piezoelectric response associated 
with such a homogeneous deformation, as already illus- 
trated via Eq. (§|). 

However, for the same situation, the improper piezo- 
electric tensor would have nonzero elements. For the cho- 
sen surface termination, the crystallite has a total dipole 
moment d = N 3r yec&, and a polarization P = d/N 3 a 2 c = 
( r ye/a 2 )z as expected. Clearly this P is invariant with re- 
spect to an elongation of the crystallite along the z axis 
(strain component e zz ), but not to an elongation along 
the x or y axes (e yy or e zz ), thus explaining why there 
is a correction to c zxx but not to c zzz in Eq. (|l6|). Simi- 
lar considerations applied to shear strains and rotations 
explain the remaining entries in Eq. (ttq). 



proposed. Instead of first computing the improper re- 
sponse and then the needed corrections, the proper re- 
sponse is computed directly from Eq. ( p2~| ) or (|24|). It is 
thus clarified that the central quantities needed to deter- 
mine the proper piezoelectric response are just the vari- 
ations of the Berry phases with deformation. 
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IV. DISCUSSION 



As is evident from Eqs.(E^) and (|T^), the distinction 
between the proper and improper piezoelectric tensor is 
only present if a spontaneous polarization is present. If 
the spontaneous polarization is small, as for wurtzite 
semiconductors H], it may be a good approximation to 
neglect the corrections to the improper tensor compo- 
nents. Alternatively, linear-response methods can be 
used to compute the proper piezoelectric response di- 
rectly |J. However, for a finite-difference calculation of 
the proper piezoelectric response of a ferroelectric mate- 
rial, it is essential to take the corrections to the improper 
response explicitly into account, as was done in Ref. M. 



V. SUMMARY 



In this work, a simple and straightforward method for 
computing the proper piezoelectric response has been 
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